//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
// Implementacao das funcoes da biblioteca c++ para programa de modelagem de
// curva de luz. Versao C++ da biblioteca roche.h para modelar curvas de luz
// da estrela secundaria em Variaveis Cataclismicas. Esta verso deve ser mais
// generica e permitir outras geometrias de forma mais simples que a versao
// C. Por exeplo, calcular curva de luz para diferentes geometrias, disco,
// estrela esferica, estrela preenchendo lobo de roche e estrela nao 
// preenchendo o lobo de roche. 
//
// Os cabecalhos das funcoes se encontram no arquivo 'light_curve.h'. Com
// ele deve ser possivel de entender o funcionamento da biblioteca sem a 
// nescessidade de ler os codigos deste arquivo.
//
// == Log de atividades ==
//
// 26/05/2006 
//    -> inicio da programacao 
//    -> inclusao de comentarios 
//
// Tiago Ribeiro -- LinuxC++ -- UFSC -- 2006
//+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

#include "light_curve.h"
 
using namespace std;

Light_Curve::Light_Curve(const vector<vector<double> > &data_par)
{
  // a entrada eh um vetor de vetores... caso o modelo precise de mais de
  // um tipo diferente de parametro! por exemplo; parametros orbitais, 
  // parametros atmosfericos, parametros sei la mais o que...

  //copia o vetor de parametros para o vetor par da classe!!!
  _vorb=true;
  for(int i = 0; i < data_par.size(); i++)
    par.push_back(data_par[i]);

  //eclipse=default_ec();
}

surface Light_Curve::Star(void)
{

  return star;

}

int Light_Curve::size(void)
{
  return star.theta.size();
}

double Light_Curve::default_ec(int pos,const vector<double> &par)
{
  return 1.0;
}

double Light_Curve::flux(const double fase, bool eclipsa)
{
  double n_s, n_s2, u1, lt_o, rlob, x1, y1, z1, z2, xp, yp, zp, f2, f1;
  double x, y, z, t;
  int k, nobs = int(par[1][2]);

  k = 0;
  lt_o = 0.;
  //u1 = fase*2.*M_PI;
  u1 = fase*2.*M_PI;
  
  z1 = cos(par[0][0]);
  z2 = sin(par[0][0]);
  f1 = sin(u1);
  f2 = cos(u1);
  
  x1 = z2*f2; //sin(par[0][0])*cos(u1);
  y1 = f1*z2; //sin(u1)*sin(par[0][0]);
 
  //printf("%f", rl1-lagrangeana(1/d[1]));

  for (int i = 0; i < star.theta.size(); i++){
    
    n_s = star.nx[i]*x1+star.ny[i]*y1+star.nz[i]*z1;
    //    n_s2 = star.px[i]*x1+star.py[i]*y1+star.pz[i]*z1;

    n_s = n_s/sqrt(star.nx[i]*star.nx[i]+star.ny[i]*star.ny[i]+star.nz[i]*star.nz[i]);
 
    if(n_s > 0){ 
      lt_o += star.I[i]*n_s*obscurecimento_de_borda(n_s,nobs);
    }
  }
  // Aplicando fator de escala
  
  if(par[1][0] > 0) lt_o *= par[1][0];  
 
  //  fazer isso em outro lugar!!!
  // lt_o = lt_o+d[4]; // d[4] = brilho do disco..
  
  return lt_o;
  
}





/*
double Light_Curve::flux(const double fase,const double *(eclipsa)(const vector<double> &,const vector<double> &))
{
  double n_s, n_s2, u1, lt_o, rlob, x1, y1, z1, z2, xp, yp, zp, f2, f1;
  double x, y, z, t;
  int k, nobs = int(par[1][1]);

  k = 0;
  lt_o = 0.;
  u1 = fase*2.*M_PI;
  x1 = sin(par[0][0])*cos(u1);
  y1 = sin(u1)*sin(par[0][0]);
  z1 = cos(par[0][0]);
  z2 = sin(par[0][0]);
  f1 = sin(u1);
  f2 = cos(u1);
 
  //printf("%f", rl1-lagrangeana(1/d[1]));

  for (int i = 0; i < star.theta.size(); i++){
    
    n_s = star.nx[i]*x1+star.ny[i]*y1+star.nz[i]*z1;
    n_s2 = star.px[i]*x1+star.py[i]*y1+star.pz[i]*z1;

    n_s = n_s/sqrt(star.nx[i]*star.nx[i]+star.ny[i]*star.ny[i]+star.nz[i]*star.nz[i]);
 
    if(n_s > 0 && n_s2 > 0) {
      vector<double> vn1 (3);
      vn1[0] = x1;
      vn1[1] = y1;
      vn1[2] = z1;      
      vector<double> vp0 (3);
      vp0[0] = star.nx[i];
      vp0[1] = star.ny[i];
      vp0[2] = star.nz[i];
      lt_o += star.I[i]*n_s*obscurecimento_de_borda(n_s,nobs)* *eclipsa(vn1,vp0);
    }
  }

  // Aplicando fator de escala
  
  if(par[1][0] > 0) lt_o *= par[1][0]; 
 
  //  fazer isso em outro lugar!!!
  // lt_o = lt_o+d[4]; // d[4] = brilho do disco..
  
  return lt_o;
  
}
*/

Histograma Light_Curve::rad_prof(const double fase, const vector<double> &format)
{
  double n_s, u1, rlob, x1, y1, z1, z2, xp, yp, zp, f2, f1;
  double x, y, z, t, w, a, r2;
  int k;
  vector<vector<double> > rad_p(2);

  using namespace Fis_const;

  k = 0;
  //u1 = fase*2.*M_PI;
  u1 = fase*2.*M_PI;
  x1 = sin(par[0][0])*cos(u1);
  y1 = sin(u1)*sin(par[0][0]);
  z1 = cos(par[0][0]);
  z2 = sin(par[0][0]);
  f1 = sin(u1);
  f2 = cos(u1);

  w = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4*M_PI*M_PI*w*w), 1./3.);

  //r2 = dd[2]/(dd[1]+dd[2]);//raio orbital da sec. em unidade de a
  r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a


  //printf("%f", rl1-lagrangeana(1/d[1]));

  for (int i = 0; i < star.theta.size(); i++){
    n_s = star.nx[i]*x1+star.ny[i]*y1+star.nz[i]*z1;

    n_s = n_s/sqrt(star.nx[i]*star.nx[i]+star.ny[i]*star.ny[i]+star.nz[i]*star.nz[i]);

    if(n_s > 0){
      // guarda velocidade e Is de cada pixel visivel em rad_p
      rad_p[0].push_back((x1*star.vx[i])+y1*(star.vy[i]));
      rad_p[1].push_back(star.I[i]*n_s*obscurecimento_de_borda(n_s,int(par[1][2])));
 
      //cout << (x1*star.vx[i])+(y1*star.vy[i])-2.*M_PI*w*a*r2*f1 << " "
      //<< star.I[i]*n_s*obscurecimento_de_borda(n_s,0) << endl;
    }
  }
  
  //se format.size == 3 modelo padrao...
  // caso contrario format contem um grid de histograma...
  //if(format.size() == 3) {
  double plen = sqrt(pix_size()); //1e5
  Histograma rad_prof(rad_p, format, plen*2.);
  
  if(par[1][0] > 0) rad_prof*=par[1][0]; //*(rad_prof.max()-rad_prof.min())/rad_prof.step_med());
  
  //else rad_prof*=(rad_prof.max()-rad_prof.min())/rad_prof.step_med();

  return rad_prof;

  //if(par[1][0] > 0) rad_prof*=(par[1][0]*(format[1]-format[0])/format[2]);
  //else rad_prof*=((format[1]-format[0])/format[2]);


    //}
    //else{
    //Histograma rad_prof(format);
    
    
    //}


  // Aplicando fator de escala
  
  // if(par[1][0] > 0) lt_o *= par[1][0]; 
 
  //  fazer isso em outro lugar!!!
  // lt_o = lt_o+d[4]; // d[4] = brilho do disco..
  
}

//fazendo perfil de modelo igual aos dados...

Histograma Light_Curve::rad_prof(const double fase, const vector<double> &format, const vector<int> &pixmap)
{
  double n_s, u1, rlob, x1, y1, z1, z2, xp, yp, zp, f2, f1;
  double x, y, z, t, w, a, r2;
  int k;
  vector<vector<double> > rad_p(2);

  using namespace Fis_const;

  k = 0;
  //u1 = fase*2.*M_PI;
  u1 = fase*2.*M_PI;
  x1 = sin(par[0][0])*cos(u1);
  y1 = sin(u1)*sin(par[0][0]);
  z1 = cos(par[0][0]);
  z2 = sin(par[0][0]);
  f1 = sin(u1);
  f2 = cos(u1);

  w = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4*M_PI*M_PI*w*w), 1./3.);

  //r2 = dd[2]/(dd[1]+dd[2]);//raio orbital da sec. em unidade de a
  r2 = 1./(par[0][2]/par[0][1]+1.);//raio orbital da sec. em unidade de a
  double q = par[0][2]/par[0][1];
  double rl2 = 0.49*pow(q,2./3.)/(0.69*pow(q,2./3.)+log(1+pow(q, 1./3.)))*a; //6.955e10;

  // Calculo o tamanho do pixel em km/s 

	double plen = sqrt(pix_size()); //1e5
 //cout << plen/1e5 << endl;
//	exit(0);
  //printf("%f", rl1-lagrangeana(1/d[1]));
  //vector<int> lpmap;

  //lpmap.reserve(pixmap.size()); // crio novo mapeamento de pixels para contar pixels que caem fora dos seus limites

  for (int i = 0; i < star.theta.size(); i++){
    n_s = star.nx[i]*x1+star.ny[i]*y1+star.nz[i]*z1;

    n_s = n_s/sqrt(star.nx[i]*star.nx[i]+star.ny[i]*star.ny[i]+star.nz[i]*star.nz[i]);

    if(n_s > 0){
      // guarda velocidade e Is de cada pixel visivel em rad_p
      rad_p[0].push_back((x1*star.vx[i])+y1*(star.vy[i]));
      rad_p[1].push_back(star.I[i]*n_s*obscurecimento_de_borda(n_s,int(par[1][2])));
      //lpmap.push_back(pixmap[i]);
      //cout << (x1*star.vx[i])+(y1*star.vy[i])-2.*M_PI*w*a*r2*f1 << " "
      //<< star.I[i]*n_s*obscurecimento_de_borda(n_s,0) << endl;
    }
    else {
      //a funcao load_data vai cuidar de nao contar esses pixels...
      rad_p[0].push_back((x1*star.vx[i])+y1*(star.vy[i]));
      rad_p[1].push_back(0.0); //por via das duvidas I = 0;
    }
  }

  

  //carrego grid
  Histograma rad_prof(format);

  //cout << lpmap.size() << " " << pixmap.size() << endl;

  //carrego dados
	rad_prof.load_data(rad_p, pixmap,plen*2.);

  if(par[1][0] > 0) rad_prof*=par[1][0];// *(rad_prof.max()-rad_prof.min())/rad_prof.step_med());
  //else rad_prof*=(rad_prof.max()-rad_prof.min())/rad_prof.step_med();

  return rad_prof;
  
}


//devolve velocidade do pixel pix

double Light_Curve::PixVel(int pix, double fase)
{
  double u1 = fase*2.*M_PI;
  double x1 = sin(par[0][0])*cos(u1);
  double y1 = sin(u1)*sin(par[0][0]);
  double z1 = cos(par[0][0]);
  double z2 = sin(par[0][0]);
  double f1 = sin(u1);
  double f2 = cos(u1);

  using namespace Fis_const;

  double w = 1./(par[0][3]*3600.); //frequencia angular
  //sep. orbital
  double a = pow(G_CGS*MSOL_CGS*(par[0][1]+par[0][2])/(4*M_PI*M_PI*w*w), 1./3.);
  //raio orbital da sec. em unidade de a
  double r2 = 1./(par[0][2]/par[0][1]+1.);

  return (x1*star.vx[pix])+(y1*star.vy[pix]);
}

int Light_Curve::PixVelPos(double vel, const vector<double> &format)
{
  if(format.size() < 3) {
    cout <<"Tamanho de vetor nao pode ser menor que 3!\n";
    exit(0);
  }
  int pixpos = 0;

  if(vel < format[0]) return pixpos;

  for(double dvel = format[0]; dvel < format[1]; dvel+=format[2]){
    //cout << dvel-format[2]/2 << " " 
    // << dvel+format[2]/2 << " " 
    // << vel << " " 
    // << pixpos << endl;
    if(((dvel-format[2]/2. <= vel) && (vel < dvel+format[2]/2.))){
      return pixpos;
    }
    pixpos++;
  }

  cout << "# - Vmin : " << format[0]-format[2] << endl
       << "# - Vmax : " << format[1]-format[2] << endl
       << "# - Velocidade: " << vel << endl;
  exit(0);
}

double Light_Curve::obscurecimento_de_borda(const double nv, const int law)
{
  //seleciona lei de obscurecimento de borda!
  // quando incluir uma lei lembrar de incluir na definicao de numero de 
  // parametros. int Light_Curve::npar_obs_borda(const int law) mais abaixo.

  switch(law) {
  case -1:
    return 1;
  case 0:
    return sqrt_law(nv); // lei da raiz quadrada
  case 1:
    return lin_law(nv); // Lei Linear
  case 2:
    return quad_law(nv); // Lei quadratica
  case 3:
    return ln_law(nv); // Lei logaritmica
  case 4:
    return nlin_claret_law(nv); // Lei proposta por Claret A. 2000AA.363..1081
  default:
    cout << "Leis de obscurecimento de borda definidas:\n"
      	 << "\t-1 - nenhuma lei (retorna 1).\n"
	 << "\t0 - lei da raiz quadrada.\n"
      	 << "\t1 - lei linear.\n"
      	 << "\t2 - lei quadratica.\n"
	 << "\t3 - lei Logaritmica.\n"
	 << "\t4 - Lei proposta por Claret A. 2000AA.363..1081\n";
    exit(1);
  }

}

int Light_Curve::npar_obs_borda(const int law)
{
  //devolve numero de parametros para lei de obscurecimento de borda!

  switch(law) {
  case -1:
    return 0;
  case 0:
    return 2; // lei da raiz quadrada
  case 1:
    return 1; // Lei Linear
  case 2:
    return 2; // Lei quadratica
  case 3:
    return 2; // Lei logaritmica
  case 4:
    return 4; // Lei proposta por Claret A. 2000AA.363..1081
  default:
    cout << "Leis de obscurecimento de borda definidas:\n"
	 << "\t-1 - nenhuma lei (retorna 1).\n"
	 << "\t0 - lei da raiz quadrada.\n"
	 << "\t1 - lei linear.\n"
      	 << "\t2 - lei quadratica.\n"
	 << "\t3 - lei Logaritmica.\n"
	 << "\t4 - Lei proposta por Claret A. 2000AA.363..1081\n";
    exit(1);
  }

}

double Light_Curve::sqrt_law(const double nv)
{
  const int N_OBS_PAR = 2;
  
  if(par[2].size() != N_OBS_PAR){
    cerr << "ERRO: Numero de parametros para lei de obscurecimento de borda\n"
	 << "ERRO: errado. Terminando execucao do programa!\n";
    exit(1);
  }

  return (1-par[2][0]*(1-nv)-par[2][1]*(1-sqrt(nv)));
}

double Light_Curve::lin_law(const double nv)
{
  const int N_OBS_PAR = 1;
  
  if(par[2].size() != N_OBS_PAR){
    cerr << "ERRO: Numero de parametros para lei de obscurecimento de borda\n"
	 << "ERRO: errado. Terminando execucao do programa!\n";
    exit(1);
  }

  return (1-par[2][0]*(1-nv));
}

double Light_Curve::quad_law(const double nv)
{
  const int N_OBS_PAR = 2;
  
  if(par[2].size() != N_OBS_PAR){
    cerr << "ERRO: Numero de parametros para lei de obscurecimento de borda\n"
	 << "ERRO: errado. Terminando execucao do programa!\n";
    exit(1);
  }

  return (1-par[2][0]*(1-nv)-par[2][1]*(1-nv)*(1-nv));
}

double Light_Curve::ln_law(const double nv) 
{
  const int N_OBS_PAR = 2;
  
  if(par[2].size() != N_OBS_PAR){
    cerr << "ERRO: Numero de parametros para lei de obscurecimento de borda\n"
	 << "ERRO: errado. Terminando execucao do programa!\n";
    exit(1);
  }

  return (1-par[2][0]*(1-nv)-par[2][1]*nv*log(nv));
}

double Light_Curve::nlin_claret_law(const double nv)
{
  const int N_OBS_PAR = 4;
  
  if(par[2].size() != N_OBS_PAR){
    cerr << "ERRO: Numero de parametros para lei de obscurecimento de borda\n"
	 << "ERRO: errado. Terminando execucao do programa!\n";
    exit(1);
  }

  double i_nv = 1;

  for(int i = 0; i < N_OBS_PAR; i++) 
    i_nv -= par[2][i]*(1-pow(nv,double(i/2.)));

  return i_nv;
}

void Light_Curve::spot_maker_test(void)
{
	if(star.I.size() != star.theta.size()) {
		cerr << "WARNING: Intensidades ainda nao definidas!\n"
		<< "WARNING: SPOT_MAKER nao pode continuar!\n"
		<< "WARNING: Voltando para programa...\n";
		return void();
	}	

	//cout << par[3][0] << " " << par[3][1] << endl;
	for(int i = 0; i < star.theta.size(); i++)
	{
		for(int ns = 0; ns < int(par[1][3]); ns++){
			
			bool test_theta = ( ( star.theta[i] > (par[3+ns][0]-par[3+ns][2])*M_PI/180.0 ) && (star.theta[i] < (par[3+ns][0]+par[3+ns][2])*M_PI/180.0  ) );
			test_theta = test_theta || ( ( star.theta[i] - 2.*M_PI > (par[3+ns][0]-par[3+ns][2])*M_PI/180.0 ) && (star.theta[i] - 2.*M_PI < (par[3+ns][0]+par[3+ns][2])*M_PI/180.0  ) );

			bool test_phi = ( ( star.phi[i] > (par[3+ns][1]-par[3+ns][3])*M_PI/180. ) && (star.phi[i] < (par[3+ns][1]+par[3+ns][3])*M_PI/180.) ) ;
			test_phi = test_phi || ( ( star.phi[i] + M_PI > (par[3+ns][1]-par[3+ns][3])*M_PI/180. ) && (star.phi[i] + M_PI < (par[3+ns][1]+par[3+ns][3])*M_PI/180.) ) ;
			
			if ( test_theta && test_phi ) {
				sum_flux_pix(i, par[3+ns][4]);
			}
		}
		// Para pixel, calcula distancia ate centro do spot
		//for(int ns = 0; ns < int(par[1][3]); ns++){
			//double cosDs = cos(M_PI/2.-star.theta[i])*cos(M_PI/2.-par[3+ns][0]*M_PI/180.0) + sin(M_PI/2.-star.theta[i])*sin(M_PI/2.-par[3+ns][0]*M_PI/180.0)*cos(star.phi[i] - par[3+ns][1]*M_PI/180.0);
			//sum_flux_pix(proc_pix(par[3+ns][0]*M_PI/180.0, par[3+ns][1]*M_PI/180.0), par[3+ns][4]);
			//if (acos(cosDs) < par[3+ns][2]*M_PI/180.0) {
			//sum_flux_pix(i, par[3+ns][4]*exp(-(cosDs*cosDs)/(par[3+ns][3]*M_PI/180.)));
			//	sum_flux_pix(i, par[3+ns][4]);
			//}
			//	}
	}

}

void Light_Curve::spot_maker(void)
{
  // colocar spots na superficie da estrela

  if(star.I.size() != star.theta.size()) {
    cerr << "WARNING: Intensidades ainda nao definidas!\n"
	 << "WARNING: SPOT_MAKER nao pode continuar!\n"
	 << "WARNING: Voltando para programa...\n";
    return void();
  }
  /*
  const double D2P = M_PI/180;

  for(int i = 0; i < star.theta.size(); i++){

    for(int ns = 0; ns < int(par[1][3]); ns++){
      if ( ((star.theta[i] >= D2P*(par[3+ns][1]-par[3+ns][3])) && (star.theta[i] <= D2P*(par[3+ns][1]+par[3+ns][3])) && (star.phi[i] >= D2P*(par[3+ns][1]-par[3+ns][3])) && (star.phi[i] <= D2P*(par[3+ns][1]+par[3+ns][3]))) ){
	sum_flux_pix(i,par[3+ns][4]);
	//cout << "Fazendo spot...\n";
      }
    }
  }
  
  for(int ns = 0; ns < int(par[1][3]); ns++)
    cout << par[3+ns][1] << " " << D2P*par[3+ns][3] 
	 << " " << par[3+ns][4] << endl;

 ----------------------------------
    double cur_pp = star.theta[i];
    double cur_pt = star.phi[i];
    double cur_I = 0.;
    for(int ns = 0; ns < int(par[1][3]); ns++){
      double pu = 1; //(cur_pt-(par[3+ns][1]*M_PI/180.)-M_PI/2.)/(par[3+ns][3]*M_PI/180.);
      double pb = sin(par[3+ns][1]*M_PI/180)*(asin(sin(cur_pp))-(par[3+ns][1]*M_PI/180.))/(par[3+ns][3]*M_PI/180.);
      cur_I += par[3+ns][4]*exp(-(pu*pu+pb*pb));
      //cout << cur_pt << " " << par[3+ns][1] << " " << par[3+ns][3] << endl;
    }


  */

  /*
  if(int(par[1][3]) == 1)
    cout << "#Fazendo " << int(par[1][3]) << " spot...\n";
  else
    cout << "#Fazendo " << int(par[1][3]) << " spots...\n";
  */



  for(int i = 0; i < star.theta.size(); i++){
    double cur_I = 0.;
    double px = sin(star.phi[i]);
    double py = cos(star.phi[i])*cos(star.theta[i]);
    double pz = cos(star.phi[i])*sin(star.theta[i]);
	  
    for(int ns = 0; ns < int(par[1][3]); ns++){
      //double rx = star.px[i]-sin((90-par[3+ns][1])*M_PI/180.);
      //double ry = star.py[i]-cos((90-par[3+ns][1])*M_PI/180.);
      //double rz = star.pz[i];
      
      double x0 = sin((90-par[3+ns][0])*M_PI/180.);
      double y0 = cos((90-par[3+ns][0])*M_PI/180.)*cos((par[3+ns][1])*M_PI/180.);
      double z0 = cos((90-par[3+ns][0])*M_PI/180.)*sin((par[3+ns][1])*M_PI/180.);
      
      double x1 = sin((90-par[3+ns][1]+par[3+ns][2])*M_PI/180.);
      double y1 = cos((90-par[3+ns][1]+par[3+ns][2])*M_PI/180.)*cos((par[3+ns][3])*M_PI/180.);
      double z1 = cos((90-par[3+ns][1]+par[3+ns][2])*M_PI/180.)*sin((par[3+ns][3])*M_PI/180.);
      
      double dx = px-x0;
      double dy = py-y0;
      double dz = pz-z0;
      
      double dr = sqrt(dx*dx+dy*dy+dz*dz);
      double dmax = sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)+(z0-z1)*(z0-z1));
      
		
      //cout << rx << " " << ry << " " << rz << endl;
      //if(rr < sin(par[3+ns][3]*M_PI/180))
      //if(sin(par[3+ns][3]*M_PI/180) > 0)
		if( dr < dmax){
			cur_I += par[3+ns][4]*exp(-(dr*dr)/(par[3+ns][3]/180));
		} 
		
	//cur_I += par[3+ns][4]*exp(-(dr*dr)/(par[3+ns][3]/180));
    }
    
    sum_flux_pix(i, cur_I);
  }


}

void Light_Curve::albedo(pixel& info, int mult=0)
{
  // Esta funcao calcula a influencia de um pixel com intensidade I colocado na posicao x,y,z do espaco
  // sobre a distribuicao de brilho do objeto. Esta funcao deve ser desenhada de modo a comportar o
  // calculo de multiplas iteracoes. Ou seja aquecimento da fonte pela estrela.
  // Parametros de entrada. 
  // pixel -> veja light_curve.h para definicao
  // mult -> Flag para multiplas refleccoes
  //par[1][6] = albedo
  
  for(int i = 0; i < star.theta.size(); i++){
    //1o - Progetar vetor na direcao do pixel da estrela
    double npsx= star.px[i]*info.px;
    double npsy= star.py[i]*info.py;
    double npsz= star.pz[i]*info.pz;
    double n_s=(star.nx[i]*npsx+star.ny[i]*npsy+star.nz[i]*npsz); 
    double nn_s=star.nx[i]*info.px+star.ny[i]*info.py+star.nz[i]*info.pz; 
    if(nn_s>0.0){
      double rxx = info.px-star.px[i];
      double ryy = info.py-star.py[i];
      double rzz = info.pz-star.pz[i];
      sum_flux_pix(i,n_s*info.I*par[1][6]/sqrt(rxx*rxx+ryy*ryy+rzz*rzz));
    }
  }

}

//soma f ao pixel n
void Light_Curve::sum_flux_pix(const int n,const double f)
{
  if(n < 0) {
    cout << "***WARNING: n < 0 -> Nada para fazer...\n";
    return void();
  }
  if(n >= star.I.size()){
    cout << "***WARNING: n > " << star.I.size() << " = pix_max.\n";
    return void();
  }

  // se ta tudo nos conformes soma f ao valor de I

  star.I[n]+=f;
}

//soma f ao pixel n
void Light_Curve::set_flux_pix(const int n,const double f)
{
	if(n < 0) {
		cout << "***WARNING: n < 0 -> Nada para fazer...\n";
		return void();
	}
	if(n >= star.I.size()){
		cout << "***WARNING: n > " << star.I.size() << " = pix_max.\n";
		return void();
	}
	
	// se ta tudo nos conformes soma f ao valor de I
	
	star.I[n]=f;
}

//devolve posicao de pixel mais proximo
int Light_Curve::proc_pix(const double theta, const double phi)
{

  if(star.theta.size() == 0 || star.phi.size() == 0){
    cout << "***WARNING: Superficie ainda nao definida...\n\t"
	 << "Nada para fazer.\n";
    return -1;
  }
  if(star.theta.size() != star.phi.size()){
    cout << "***ERRO: Vetores de theta e phi com tamanhos diferentes...\n\t"
	 << "fechando programa...\n";
    exit(1);
  }

  //procura pixel mais proximo!

  double df = 2*M_PI; //distancia em f
  double dp = 2*M_PI; // distancia em p
  int ppix = 0;

  for(int i = 0; i < star.theta.size(); i++){
    double cf = star.theta[i]; //current f
    double cp = star.phi[i];  //current p
    double dd = sqrt(cos(df)*cos(df)+cos(dp)*cos(dp));
    double dc = sqrt(cos(cf-theta)*cos(cf-theta)+cos(cp-phi)*cos(cp-phi));

    if(dc < dd){ // se distancia atual menor que anterior
      df = acos(cos(cf-theta)); //distancia em f igual a atual
      dp = acos(cos(cp-phi)); //distancia em p igual a atual
      ppix = i; //esse eh o pixel mais proximo...
    }
  }

  return ppix;

}

//devolve posicao de pixel mais proximo
int Light_Curve::proc_pix(const double x, const double y, const double z)
{

  if(star.px.size() == 0 || star.theta.size() == 0){
    cout << "***WARNING: Superficie ainda nao definida...\n\t"
	 << "Nada para fazer.\n";
    return -1;
  }
  if(star.px.size() != star.py.size()){
    cout << "***ERRO: Vetores de theta e phi com tamanhos diferentes...\n\t"
	 << "fechando programa...\n";
    exit(1);
  }
  if(star.px.size() != star.pz.size()){
    cout << "***ERRO: Vetores de theta e phi com tamanhos diferentes...\n\t"
	 << "fechando programa...\n";
    exit(1);
  }

  //procura pixel mais proximo!

  double dx = 2*M_PI; //distancia em x
  double dy = dx; // distancia em y
  double dz = dx; // z
  int ppix = 0;

  for(int i = 0; i < star.px.size(); i++){
    double cx = star.px[i]; //current x
    double cy = star.py[i];  //current y
    double cz = star.pz[i];  //current z

    double dd = sqrt(dx*dx+dy*dy+dz*dz);
    double dc = sqrt((cx-x)*(cx-x)+(cy-y)*(cy-y)+(cz-z)*(cz-z));

    if(dc < dd){ // se distancia atual menor que anterior
      dx = (cx-x); //distancia em x igual a atual
      dy = (cy-y); //distancia em y igual a atual
      dz = (cz-z); // z
      ppix = i; //esse eh o pixel mais proximo...
    }
  }

  return ppix;

}

//devolve numero de pixels com distancia menor do que Dist e coloca os indices dos pixels na variavel, pmap 
int Light_Curve::close_pix(const double px, const double py, const double pz, const double Dist, double*& pmap)
{
	
	if(star.px.size() == 0 || star.theta.size() == 0){
		cout << "***WARNING: Superficie ainda nao definida...\n\t"
		<< "Nada para fazer.\n";
		return -1;
	}
	if(star.px.size() != star.py.size()){
		cout << "***ERRO: Vetores de theta e phi com tamanhos diferentes...\n\t"
		<< "fechando programa...\n";
		exit(1);
	}
	if(star.px.size() != star.pz.size()){
		cout << "***ERRO: Vetores de theta e phi com tamanhos diferentes...\n\t"
		<< "fechando programa...\n";
		exit(1);
	}
	
	//procura pixel em raio R <= Dist!
	vector<double> vec_pmap(0);

	int npix = 0;
	for(int i = 0; i < star.px.size(); i++){
		double cx = star.px[i]; //current x
		double cy = star.py[i];  //current y
		double cz = star.pz[i];  //current z
		
		double dc = sqrt((cx-px)*(cx-px)+(cy-py)*(cy-py)+(cz-pz)*(cz-pz));
		
		if(dc <= Dist){ // se distancia atual menor Dist
			npix++;
			vec_pmap.push_back(i);
		}
	}
	
	pmap = new double [npix];
	
	for (int i = 0; i < npix;i++)
		pmap[i] = vec_pmap[i];
	
	return npix;
	
}


//Funcao de visibilidade

double Light_Curve::VFunc(int n, double fase)
{
  double u1 = fase*2.*M_PI;
  double x1 = sin(par[0][0])*cos(u1);
  double y1 = sin(u1)*sin(par[0][0]);
  double z1 = cos(par[0][0]);
  double n_s = star.nx[n]*x1+star.ny[n]*y1+star.nz[n]*z1;

  return n_s/sqrt(star.nx[n]*star.nx[n]+star.ny[n]*star.ny[n]+star.nz[n]*star.nz[n]);
}

